Towards Safe and Effective High-Order Legendre Transforms with Applications to FFTs for the 2-sphere

This paper proposes and investigates a method for reducing potential numerical instability in the fast spherical harmonic transform algorithms proposed previously in [5,11]. The key objective of this study is a numerically reliable fast algorithm for computing the discrete Legendre transform (DLT); that is, the projection of sampled data onto the associated Legendre functions within a specified range of degrees. A simple divide-and-conquer approach derives from a factorization of high-degree Legendre functions into Legendre functions of lower-degree, exploiting the fact that the complexity of projection onto Legendre functions decreases with decreasing degree. Combining the resulting fast DLT algorithms for each relevant order of associated Legendre function results in an O(N log2 N) algorithm for computing the spherical harmonic expansion of a function sampled at N points on the sphere. While fast DLT algorithms of this form are exact in exact arithmetic, actual (finite precision) implementations of the earlier variants display instabilities which generally grow with the the order of the associated Legendre function in the transform. Here we return to the basic algorithm, present a slight modification of the general schema and examine the error mechanism for the higher-order cases. This study suggests a new approach to the high-order DLTs achieved by substitution of a simple alternative factorization of the associated Legendre functions in place of that used previously. This technique improves stability significantly for a wide range of useful problem sizes and may be used with any of the variants of the basic algorithm previously proposed. We present a description of the use of the new Legendre decomposition in the basic fast algorithm along with numerical experiments demonstrating the large advances in stability and efficiency of the new approach over the previous results.